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ABSTRACT 


Specialized  nonlinear  crack  tip  finite  elements 
which  include  the  amplitude  of  the  plastic  singular  solu- 
tion as  an  additional  unknown  are  investigated  to  deter- 
mine their  capability  to  directly  predict  the  J-integral 
for  cracked  elastoplastic  bodies  without  recourse  to 
numerical  evaluation  of  the  J-integral  over  arbitrarily 
chosen  paths.  The  special  elements  are  used  in  conjunc- 
tion with  conventional  12-node  quadrilateral  isoparametric 
elements.  Power  hardening  and  multilinear  representations 
of  the  nonlinear  stress-strain  curve  are  considered  for 
both  deformation  and  incremental  theories  of  plasticity. 
Numerical  results  are  presented  which  demonstrate  that 
the  elements  fail  to  provide  accurate  direct  calculation 
of  J,  but  that  they  lead  to  improved  estimates  of  J based 
on  path  calculations.  It  is  concluded  that  special  ele- 
ments at  the  crack  tip  improve  the  accuracy  of  path  values 
of  J,  but  that  the  special  elements  themselves  predict  ac- 
curate values  of  J only  in  materials  with  high  strain 
hardening  slope. 


ADMINISTRATIVE  INFORMATION 

This  work  was  authorized  and  funded  within  the  Submarine  Structures 
Exploratory  Development  Program  under  Program  Element  62543N,  Project 
SF  43.422.592,  and  Work  Unit  1720-592. 

INTRODUCTION 

The  finite  element  calculation  of  linear  elastic  stress  intensity 

factors  for  planar  or  axisymmetric  bodies  has  become  commonplace  and  highly 

accurate.  Indirect  methods,  in  which  conventional  elements  are  forced  to 

display  the  correct  near  tip  displacement  field,  have  been  discussed  by 
1*  2 

Henshell  and  Shaw,  and  by  Barsoum  for  the  8-node  quadrilateral  isopara- 

3 

metric  element,  and  by  Pu  et  al.  for  the  12-node  quadrilateral  isopara- 
metric element.  In  these  methods,  a corner  node  corresponds  to  the  crack 
tip,  and  edge  nodes  adjacent  to  the  tip  are  moved  to  the  1/4  position  for 
the  8-node  element,  or  to  the  1/9-4/9  positions  for  the  12-node  element, 
thus  imposing  the  proper  /r  variation  of  displacement  with  respect  to 


j 

4 


*A  complete  listing  of  references  is  given  on  page  27. 
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distance  from  the  crack  tip.  Stress  intensity  factors  and  K^.  are  then 
calculated  based  upon  examination  of  nodal  displacements  in  the  vicinity 
of  the  crack  tip. 

4-6 

On  the  other  hand,  Hilton  and  Gifford  have  found  that  highly  accu- 
rate elastic  stress  intensity  factors  can  be  calculated  using  specialized 
crack  tip  elements  in  which  the  stress  intensity  factors  themselves  are 
carried  as  unknowns.  Special,  small,  circular  "core"  crack  tip  elements 
joined  along  their  periphery  to  conventional  elements  have  produced  good 
results.  An  "enrichment"  of  conventional  12-node  isoparametric  elements, 
as  suggested  originally  by  Benzley,^  has  proved  even  more  convenient  and 

has  improved  accuracy.  Both  of  these  approaches  have  been  implemented  into 

6 8 

the  APES  computer  program  * which  employs  the  12-node  isoparametric 

quadrilateral  as  a conventional  element. 

The  extension  of  embedded  singularity  finite  element  techniques  to 

the  consideration  of  elastoplastic  crack  problems  was  first  proposed  by 

9 

Hilton  and  Hutchinson.  In  this  work,  the  core  crack  tip  element  was 

employed  to  embed  the  Hutchinson^  and  Rice  and  Rosengren1'*'  (HRR)  plastic 

singular  solution  for  Mode  I crack  problems  into  a deformation  theory 

finite  element  formulation.  Using  this  approach,  the  plastic  intensity 

12 

factor  was  calculated  directly  and  the  J-integral  was  determined  from  it. 
Results  were  presented  for  the  plane  stress  problem  of  a central  crack  in 
an  infinite  domain  subject  to  a uniform  remote  tensile  stress  field  acting 
normal  to  the  crack  direction.  The  authors  extended  this  approach  to 
plane  strain  problems  and  to  the  use  of  high  order  isoparametric  elements 
instead  of  previously  used  constant  stress  triangles  remote  from  the  crack 
tip.  They  report  results  for  a number  of  geometries  and  characteristics, 
including  boundary  influence,  for  example. 

In  this  approach,  the  core  element  was  thought  to  be  particularly 
attractive  for  the  following  reasons: 

1.  It  was  not  necessary  to  select  a path,  a number  of  paths,  or  a 
"best"  path  over  which  to  evaluate  the  J-integral. 

2.  The  J-integral  was  a direct  consequence  of  the  calculation. 
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3.  The  J-integral  was  a first  order  quantity  whose  accuracy  should 
be  on  the  order  of  that  of  nodal  displacements;  evaluation  of  J over  some 
path,  on  the  other  hand,  involves  the  second  order  quantities  of  strain 
and  stress. 

Success  in  elastic  calculations  with  the  enriched  12-node  isopara- 
metric element^  was  even  more  dramatic  than  with  the  small  core  element. 

Not  only  were  problems  considerably  easier  to  idealize,  but  accuracy  of 
the  elastic  stress  intensity  factors  was  found  to  be  practically  insensi- 
tive to  the  size  of  the  enriched  elements  surrounding  the  crack  tip.  (Such 

was  not  the  case  with  the  core  element.)  This  fact  quickly  led  the 
authors  to  almost  exclusive  use  of  enriched  12-node  elements  for  elastic 
fracture  computations. 

i 

This  success,  along  with  dissatisfaction  over  the  limitations  im- 
posed by  deformatio/t  theory  plasticity  coupled  with  a power  hardening 
material  model,  has  led  the  authors  to  yet  another  approach  to  the  elasto- 
plastic  Mode  I crack  problem.  This  took  the  form  of  an  enriched  elasto- 
plastic  12-node  isoparametric  crack  tip  element,  coupled  with  conventional 
12-node  elements,  and  all  formulated  using  incremental  theory  plasticity 
and  a multilinear  representation  of  the  material  stress-strain  curve. 

From  analyzing  the  results  of  all  of  this  work,  the  authors  are  able 
to  compare  the  influences  of  the  following  in  elastoplastic  fracture 
analysis : 

1.  Different  models  of  material  behavior 

2.  Different  near-tip  elements  (core,  enriched,  1/9-4/9  and 
nonsingular) 

3.  Path  and  singularity  (direct  formulation,  first  order)  predic- 
tions of  the  J-integral 

In  the  following,  the  development  of  each  of  the  nonlinear  formulations 
will  be  outlined,  and  then  numerical  results  will  be  compared.  It  will 
be  concluded  that  while  path  calculation  of  the  J-integral  is  satisfactory 
using  singular  elastoplastic  elements,  the  directly  calculated  singular 
values  themselves  are  not,  in  general,  accurate.  Thus  the  advantage  en- 
joyed in  the  use  of  such  elements  in  the  elastic  case  is  lost  in  the 
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nonlinear  case.  The  unsatisfactory  prediction  of  singular  J values  will 
then  be  explained  and  discussion  regarding  the  need  for  special  elements 
for  accurate  path  prediction  will  be  given. 

MATERIAL  MODELS  AND  ASSOCIATED  CRACK  TIP  SINGULARITIES 
POWER  HARDENING  MATERIAL 

For  a power  hardening  material,  the  uniaxial  stress-strain  law  is 
given  by 

a for  o < a 

— y p 

o - aa  + a a(a/a  )n  1 for  a > a ^ 

yp  yp  yp 

where  a and  £ = uniaxial  stress  and  strain 
o = the  yield  stress 

yp 

E = Young's  modulus 

ot  and  n = constants  chosen  such  that  Equation  (1)  models  the  experi- 
mental stress-strain  curve 

For  such  a material  model  under  conditions  of  plane  stress,  plane  strain, 

9 

or  axial  symmetry,  Hutchinson  has  shown  that  the  near  crack  tip  fields 
are  given  for  the  Mode  I case  by 


■ V"1/n+1 5ij<e) 


(2a) 


p p -n/n+1  ~ 

Elj=Elj'“E  ' V6) 


(2b) 


p 1/n+l  ~ ,Q. 
Ui  = Uio  + a E r Ui(6) 


(2c) 
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where  o_  , e , and  = stress,  strain,  and  displacement  components 

r and  0 = polar  coordinates  centered  at  the  crack  tip  as 

shown  in  Figure  1 

Kp  = the  plastic  intensity  factor 


Figure  1 - Polar  Coordinate  System  at  Crack  Tip  and 
Contour  for  J Path  Integral 

The  functions  a„(Q),  £^(0),  and  u^(0)  are  determined  by  the  numerical 
solution  of  a nonlinear  fourth  order  ordinary  differential  equation;  they 
are  dependent  on  the  hardening  exponent  n and  are  distinct  for  the  cases 
of  plane  stress  and  plane  strain.  The  near-tip  fields.  Equations  (2), 


r 


are  limited  to  the  region  for  which  plastic  strain  components  dominate 

the  corresponding  elastic  components.  Thus  the  physical  region  over  which 

Equation  (2)  applies  is  a function  of  applied  load  and  increases  in  size 

as  the  applied  load  increases. 

12 

Rices'  J-integral  given  in  terms  of  the  strain-energy  density  W, 


is 


0. .de.  . 
iJ  ij 


(3) 


J 


I 


du . 

Wdy  - Ti  ds 


(4) 


The  notation  of  Equation  (4)  is  shown  in  Figure  1.  It  has  been  shown  that 
along  any  path  T (Figure  1),  from  the  lower  crack  face  to  the  upper,  J is 
path  independent  for  deformation  theory  plasticity  with  no  unloading. 
Further,  its  value  is  directly  related  to  the  plastic  intensity  factor  K 
when  these  conditions  are  satisfied,  i.e.. 


J 


(5) 


where  I is  a definite  integral  dependent  only  on  the  hardening 
exponent  n. 


BILINEAR  AND  MULTILINEAR  MATERIALS 

The  associated  near-tip  field  found  by  Hutchinson  for  bilinear 

material  models  has  not  been  employed  previously  in  conjunction  with  the 

finite  element  method  to  determine  J values,  etc.  For  this  reason,  the 

description  here  will  contain  more  detail  than  that  for  the  power  harden- 

4 

ing  material  fully  described  elsewhere. 
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The  uniaxial  multilinear  model  of  material  behavior  is 


£ = f + T (ar°yp)  + T (°2"ai)  +'*‘  t (°-a 


m-r 


(6) 


where  o , < a < 0 and  a is  defined  as  (EAe  -A a )/Ao  . The  number  of 
m-1  — m m m m m 

inelastic  plastic  segments  in  the  material  model,  N,  is  greater  than  or 

equal  to  m.  The  cutoff  of  the  last  segment  is  at  infinity,  i.e. , £ = 00 . 

10  N 
The  Hutchinson  work  is  based  on  the  special  case  of  bilinear 

material  response,  N = 1,  and  the  assumption  that  the  region  immediately 

surrounding  the  crack  tip  is  yielded.  He  obtains  the  asymptotic  solution 

for  the  Airy  stress  function 


= 


/H 


(cos  0/2+1/3  cos  30/2) 


(7) 


applicable  to  both  plane  stress  and  plane  strain.  Thus  the  plastic  singu- 
lar solution  for  the  near  tip  stress  components  has  the  same  (r,0)  depend- 
ence as  in  the  elastic  case  within  the  yielded  zone. 

The  generalization  of  Equation  (6)  for  multiaxial  stress  states  is 


£.  . 

ij 


(1+v) 

E 


'ij 


V 

E 


a 6 . 

PP  ij 


V[(ai0yp)  + (Vai)0l  +- 


(a  -a  , )o 

m m-1  m- 


l]/°e 


ij 


(8) 


when  (a  ) , < a < (a  ) . In  the  above,  v is  Poisson's  ratio,  5..  is  the 

e m-1  e — e m ij 

Kroneker  delta,  Oe  is  the  effective  or  von  Mises  stress,  s is  the  devi- 
atoric  sfess  tensor,  and  the  are  the  parameters  defining  the  slopes  of 
the  various  segments  describing  the  hardening  portions  of  the  uniaxial 
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r 


stress-strain  curve.  Using  Equation  (8),  the  asymptotic  stress,  strain, 
and  displacement  fields  within  the  yielded  zone  may  be  obtained  from  the 
Airy  stress  function.  The  results  for  the  displacement  fields  are: 


for  plane  stress  and 


u = u + 


j.  fi F 

4E  V 7T 


(5-3v+  | aN-83)  cos  | - (l+\H-  | aN)  cos 


+ 0 (r ) 

(10a) 


_£ 

4E 


# [(7_vf  T V83)  sin  I - (1+XH_  I %)  sin  T 


+ 0(r) 


(10b) 


for  the  plane  strain  case  with  6 = (\H-a^/2)  /(1+a^). 

The  J-integral  can  be  expressed  in  terms  of  the  plastic  intensity 
factor  Kp  by  carrying  out  the  integration  on  a path  in  the  region  governed 
by  the  plastic  singularity  solution.  Equations  (9)  or  (10).  The  results 
are 


1+Y  /2-1.5vy+a, 


2 

J = (1+a^)  K^/e  (plane  stress) 

(i_  ii  + 1! 

N \ 4 + 2 


Kp/E  (plane  strain) 


(11a) 


(lib) 


with  Y = (v+ol /2)  • (1+a  ) 
N N 
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It  is  important  to  note  that  the  development  of  Equations  (9)  and  (10) 
differs  significantly  from  the  corresponding  derivation  of  Equations  (2) 
for  the  power  hardening  material  in  that: 

1.  Equations  (9)  and  (10)  contain  the  elastic  and  the  plastic  con- 
tributions to  the  asymptotic  displacement  field  while  Equations  (2)  contain 
only  the  plastic  contribution. 

2.  In  Equations  (9)  and  (10),  the  contributions  of  order  0(r)  which 
are  neglected  contain  the  portion  of  the  displacement  component  derived 
from  the  s^j/°e  term  in  Equation  (8).  In  other  words,  the  asymptotic  solu- 
tion has  been  found  based  on  the  strain-stress  assumption 


_ (1+v)  „ v , , 3 

Eij  E ij  E °pp6ij  + 2E  aNSij 


(12) 


SPECIALIZED  CRACK  TIP  ELEMENTS 

CORE  ELEMENT 

4 

A semicircular  core  element  has  been  developed  for  elastoplastic 
Mode  I problems  in  which  the  assumed  displacement  field  is  given  by  the 
asymptotic  solution.  Equations  (2),  (9),  or  (10),  with  unknown  parameters 
uq  (the  x displacement  component  of  the  crack  tip)  and  K^.  This  small 
special  element  is  connected  to  standard  12-node  isoparametric  elements 
along  its  periphery  by  using  the  last  of  Equations  (2),  (9),  or  (10)  as 
displacement  constraints  for  nodes  falling  on  tne  boundary.  The  element 
thus  has  a variable  number  of  nodal  points;  usually,  three  or  four  standard 
12-node  elements  adjoin  its  boundary.  The  related  finite  element  calcula- 
tions result  in  direct  predictions  of  K and  u as  well  as  all  nodal  dis- 

P o 

placement  components,  strains,  and  stresses.  Employing  Equations  (5), 
(11a),  or  (lib),  the  authors  were  able  to  directly  calculate  the  J-integxal 
as  an  unknown  associated  with  the  plastic  singular  solution.  In  addition, 
provision  has  been  made  to  evaluate  J along  a path  near  the  crack  tip  in 
the  first  ring  of  conventional  elements  around  the  core  element  according 
to  Equation  (4). 
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ENRICHED  ELEMENT 


The  concept  of  enriching  the  conventional  element  displacement  assump- 
tion with  the  asymptotic  displacement  field  appropriate  for  fracture  analy- 
sis originated  with  Benzley^  and  has  been  extensively  applied  for  the 
elastic  case  by  the  authors. The  application  described  here  for  the 
elastoplastic  Mode  1 crack  problem  is  new.  In  the  approach  the  multilinear 
material  model  is  used  and  leading  term  of  the  asymptotic  expansion  for  the 
near-tip  crack  displacement,  Equations  (9)  or  (10),  is  added  to  the  usual 
polynomial  displacement  field  within  the  element.  The  present  application 
involves  enriching  12-node  isoparametric  elements,  i.e. , 


u 


a,  + a.s  + a.t  + a.s 
12  3 4 


+ . . . 


a,  „st  + K u(s,  t) 

12  p 


or 


(13) 


u = [P] {a}  + K u(s, t) 

P 

where  the  "a's"  are  generalized  displacements,  (s,t)  are  the  local  coor- 
dinates associated  with  the  isoparametric  element,  and  u(s,t)  represents 
the  term  (transformed  to  local  element  coordinates)  in  Equations  (9)  or  (10) 
in  the  expression  for  u whose  coefficient  is  Kp.  The  expression  for  the  v 
component  of  displacement  is  similar. 

Equation  (13)  is  evaluated  at  each  of  the  twelve  nodes  to  obtain  the 
set  of  equations  for  the  nodal  values  of  u,  i.e.,  {u},  in  terms  of  the 
vector  {a}  as 


{u}  = [C] {a}  + Kp{2}  (14) 

Inversion  yields 

{a}  = [C]"1  ({u}  - Kp{Q})  (15) 
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Back  substitution  into  Equation  (13)  gives 


u = [Pile]”1  ( {u}-K  {u})  + K Q 
P P 


(16) 


But  from  this  equation,  it  is  clear  that  [P][C]  1 is  nothing  more  than  the 

set  of  usual  isoparametric  shape  functions  [N]  = N12^‘ 

Therefore 


u 


E 

^ =i 


N. ( s , t ) u , + K 
i i p 


u(s,t)- 


E 


Ni(s,t)u1 


(17) 


and,  similarly. 


v 


E 

■i=i 


N.  (s,t)v.  + K 
i i p 


v(s,t)- 


E 

-i  = 1 


Ni(s,t)vi 


(18) 


Equations  (17)  and  (18)  are  the  enriched  element  displacement  approxima- 
tions, giving  the  displacement  components  within  the  element  in  terms  of 
their  nodal  values  and  the  plastic  intensity  factor  K . The  enriched 
element  stiffness  matrix  is  developed  from  these  equations  in  the  usual 
manner,  with  the  exception  that  higher  order  Gauss  quadrature  (8x8)  is 
employed  in  the  area  integration  to  properly  handle  the  steep  gradients 
near  the  crack  tip. 

Implementation  of  enriched  elements  of  this  type  into  a standard 
finite  element  code  is  not  as  straightforward  as  for  the  core  element  dis- 
cussed previously.  The  enriched  element  stiffness  matrix  is  of  order  25 
for  the  Mode  I case,  corresponding  to  24  unknown  nodal  displacement  com- 
ponents and  the  singular  solution  coefficient  K . Allowance  in  the  con-  \ 

P 

struction  and  solution  of  the  overall  stiffness  matrix  must  be  ™*de  for 
an  additional  row  and  column  which  correspond  to  K . 
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The  J-integral  is  calculated  from  using  Equation  (11a)  or  (lib)  as 
described  previously.  In  addition,  the  J-integral  may  be  evaluated  on  up 
to  ten  different  paths  specified  by  the  analyst  according  to  Equation  (4). 

ELEMENT  WITH  ADJUSTED  NODE  POSITIONS  (1/9-4/9) 

14 

Barsoum  has  shown  that  when  midside  nodes  of  8-node  quadrilateral 

isoparametric  elements  are  moved  to  the  quarter  point  nearest  the  crack 

tip,  a l//r  singularity  in  the  strain  field  is  developed.  Following  this 

3 

work,  Pu  and  Hussain  found  that  the  same  effect  could  be  achieved  with 
12-node  quadrilateral  isoparametric  elements  by  moving  intermediate  nodes 
to  the  1/9  and  4/9  positions  nearest  the  crack  tip.  The  writers  have  per- 
formed elastic  calculations  with  1/9-4/9  quadrilateral  elements  (rather 
than  corresponding  collapsed  triangular  elements  currently  preferred  by 
most  investigators)  and  have  obtained  highly  accurate  results.  Benzley'*'"’ 
has  shown  that  the  use  of  quarter  point  8-node  elements  in  elastic-plastic 
computations  leads  to  improved  path  values  of  the  J-integral.  At  the  same 
time,  the  plastic  singular  solution  associated  with  a multilinear  material 
model  contains  the  square  root  singularity  associated  with  the  1/9-4/9 
element.  Thus  it  became  logical  to  employ  the  1/9-4/9  element  with  a 
multilinear  material  model  as  an  additional  plastic  singularity  element. 

NUMERICAL  EXAMPLES  OF  J (PATH)  VERSUS  J (SINGULARITY) 

Based  on  the  material  described  previously,  the  authors  are  presently 
able  to  make  the  following  calculations:* 

1.  Calculation  of  J from  the  plastic  singular  solution  (hereafter 
called  J^)  of  Equation  (2)  and  calculation  of  J along  a path  (hereafter 
called  Jp)  near  the  crack  tip  for  a power  hardening  material  using  de- 
formation theory  plasticity  and  a small  semicircular  core  crack  tip 
element . 

2.  The  same  as  above  but  for  a bilinear  material  model. 


*These  calculations  have  been  made  based  on  elastic-plastic  extensions 

^ g 

of" the  AFES'xomputer  program.  ’ This  nonlinear  version  of  the  program  has 
been  given  the  acronym  PAPST  for  Plastic  Axisymme trie/ Planar  Structures . 

The  nonlinear  program  is  still  under  development  and  has  not  yet  been 
formally  documented. 
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3.  Calculation  of  J and  J for  bilinear  and  multilinear  material 

s p 

models  using  incremental  plasticity  theory  and  the  enriched  element  whose 
displacement  field  is  given  by  Equations  (17)  and  (18). 

4.  Calculation  of  J as  immediately  above,  but  with  the  use  of  the 
1/9-4/9  elements  at  the  crack  tip. 

5.  Again  the  same,  but  with  no  special  treatment  of  the  crack  tip. 

The  hypothetical  edge-cracked  specimen  whose  symmetric  upper  half  is 

shown  in  Figure  2 is  used  as  an  example  to  illustrate  the  character  of 
numerical  results  to  be  expected.  The  following  material  properties  were 
assumed:  E = 30  x 10^  psi  (20.7  * 10"*  Mpa) , V = 0.3,  and  0 = 10"*  psi 

(6900  Mpa).  The  finite  element  grid  patterns  are  also  shown  in  Figure  2 
along  with  the  contours  used  to  evaluate  J^. 

Table  1 gives  results  obtained  for  the  plane  strain  case  using  the 
core  element  approach  and  the  grid  pattern  shown  in  Figure  2b.  The  core 
radius  in  this  case  is  0.01  in.  (0.25  mm)  (0.67  percent  of  crack 
length).  Some  of  the  stress-strain  curves  for  these  calculations  are  shown 
in  Figure  3. 

The  choice  of  a small  hardening  coefficient  a(0.01)  with  a hardening 
exponent  n = 1 results  in  essentially  elastic  behavior  as  expected.  The 
Jg  and  values  in  this  range  of  material  behavior  are  in  close  agreement 
with  each  other  and  also  with  the  elastic  finite  element  predictions  (J^  ^). 

These  excellent  predictions  are  in  accordance  with  the  authors'  earlier 

. 4-6 
work. 

Reduced  material  hardening  is  handled  in  the  material  model  by  in- 
creasing the  coefficient  a and/or  the  exponent  n.  The  numerical  results 

given  in  Table  1 show  increasing  disparity  between  J and  J with  de- 

s p 

creasing  strain  hardening.  Considering  a fixed  grid  pattern  and  load 
level,  the  values  first  increase  slightly  and  then  decrease  dramatical- 
ly as  n and/or  a is  increased.  On  the  other  hand,  J increases  monoton- 

P 

ically  with  decreasing  strain  hardening.  The  observed  numerical  dependence 
of  on  a and  n is  not  physically  realistic.  The  J values,  on  the  other 
hand,  are  believed  to  exhibit  the  correct  trend. 
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Figure  2a  - Geometry  and 
Boundary  Conditions 


Figure  2b  - Nine-Element  Mesh 
(One  Core  and  Eight  Regular) 
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Figure  2d  - Six-Element  Mesh 

(Two  Enriched  and  Four 
Regular) 


Figure  2c  - Thirteen-Element 
Mesh 
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Figure  2e  - Fourteen-Element 
Mesh 


(Two  Enriched  and  Twelve 
Regular) 


Figure  2 - Geometry  and  Idealizations  of  Hypothetical 
Edge-Notched  Specimen 
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TABLE  1 - PLANE  STRAIN  CALCULATIONS  FOR  POWER  HARDENING  MODEL  FOR  VARIOUS  VALUES  OF  n AND 
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The  authors  have  previously  reported  numerical  predictions  for 
4 9 13 

J = J ’ ’ using  the  circular  core  element.  Those  results  were  for 
s ° 

materials  with  sufficient  hardening  so  that  no  significant  differences  in 
Jg  and  .1^  were  observed.  It  is  in  the  present  extension  of  this  approach 
to  lower  hardening  materials  that  the  discrepancy  between  Jg  and  became 
apparent . 

The  case  for  the  enriched  12-node  element  is  considered  next  in  the 
context  of  a bilinear  material  model  (n=i).  As  shown  in  Table  2,  when  a 
is  very  small  leading  to  essentially  elastic  conditions,  the  singular  and 
path  predictions  of  J are  again  in  close  agreement  and  also  closely  agree 
with  the  linear  elastic  finite  element  results.^  Again,  this  was  expected 
based  on  previous  work. 

Some  of  the  stress-strain  curves  assumed  in  Table  2 are  shown  in 

Figure  3.  As  the  hardening  coefficient  a increases  (leading  to  reduced 

strain  hardening),  the  difference  between  J and  J increases,  with  J 

s p s 

falling  while  increases.  As  a increases  beyond  about  0.1,  the  accuracy 

of  the  J values  rapidly  decreases.  The  J values,  on  the  other  hand, 
s p 

appear  to  be  reasonably  accurate  over  the  range  of  nonlinearity  studied. 

Table  3 presents  results  obtained  with  core  element,  enriched  ele- 
ments, 1/9-4/9  elements,  and  conventional  element  treatments  of  the  crack 
tip  singularity  for  the  cases  of  plane  strain  and  plane  stress.  The  core 
element  results  are  from  three  different  grid  patterns:  9-element  ideali- 
zation with  core  radius  rQ  of  0.02  and  0.01  in.  (0.51  and  0.25  mm),  and  a 

13-element  idealization  with  r =0.01  in.  (0.25  mm).  These  grids  are 

o 

shown  in  Figure  2.  The  elastic  results  for  J(JgSy)  exhibit  a maximum 
difference  between  grid  patterns  of  about  10  percent.  This  variation  is 
reflected  in  the  plasticity  results.  The  enriched,  1/9-4/9,  and  conven- 
tional treatments  of  the  crack  tip  shown  in  Table  3 are  for  two  grid  pat- 
terns (Figure  2)  consisting  of  6 and  14  elements.  The  elastic  results  show 
a much  weaker  dependence  of  J on  grid  pattern  than  that  found  for  the  core 

element.  Both  enriched  idealizations  provide  elastic  results  which  are 

4 

more  accurate  than  any  of  the  core  element  predictions.  Thus,  the  en- 
riched element  path  values  for  J are  believed  to  be  the  more  accurate  of 

P 

those  presented  in  Table  3. 
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*From  elastic  finite  element  solution  for 


TABLE  3 - RESULTS  OF  DIFFERENT  FINITE  ELEMENT  MODELS  F )R  BILINEAR  MATERIAL  BEHAVIOR 


integral  from  elastic  finite  element  solution  for 


The  following  general  conclusions  are  drawn  from  the  results  presented 
in  Table  3: 

1.  Values  of  J are  inaccurate  for  both  the  core  element  and  enriched 

s 

elements  for  materials  which  exhibit  low  strain  hardening. 

2.  Values  of  agree  closely  for  core  element  (deformation  theory 
plasticity)  and  enriched  element  (incremental  theory  plasticity)  idealiza- 
tions for  the  plane  stress  case,  but  disagree  moderately  in  the  plane 
strain  case. 

3.  For  plane  strain  and  plane  stress  calculations,  enriched  and 
1/9-4/9  elements  provide  about  the  same  degree  of  accuracy  for  values. 

4.  Inferior  values  of  J are  obtained  when  the 'crack  tip  is  treated 

P 

in  a nonsingular  manner.  This  trend  may  become  much  more  significant  if 
the  path  used  to  calculate  is  moved  closer  to  the  crack  tip. 

Altogether,  the  results  presented  in  Tables  1,  2,  and  3 lead  to  two 
basic  conclusions  of  fundamental  importance: 

1.  Path  values  may  be  evaluated  with  apparently  reasonable  accu- 
racy for  a wide  class  of  stress-strain  curves  using  either  the  enriched  or 
1/9-4/9  crack  tip  elements  in  conjunction  with  the  incremental  theory  of 
plasticity,  provided  that  the  chosen  path  is  not  very  close  to  the  crack 
tip.  This  conclusion  has  limited  experimental  verification*  and  has  been 
confirmed  analytically  by  comparison  with  results  obtained  in  the  ASTM- 
sponsored  round-robin  J-integral  calculation.^ 

2.  Singular  values  are  seriously  in  error  for  all  cases  except 
those  involving  materials  which  are  nearly  elastic.  This  error  occurs 
with  both  the  core  element  using  deformation  theory  plasticity  and  with 
enriched  elements  using  incremental  theory  plasticity. 

The  second  conclusion  is  in  contrast  to  the  successful  elastic  results 
which  prompted  this  work.  The  notion  that  (analogous  to  the  elastic  case) 

J could  be  calculated  directly  as  a first  order  quantity,  independent  of 
path,  and  with  accuracy  on  the  order  of  that  of  nodal  displacements,  has 
been  found  to  be  false.  The  explanation  follows. 


*This  work  has  been  reported  informally  in  DTNSRDC  Structures  Depart- 
ment Technical  Note  m-4,  "J-integral  Analysis  of  a Compact  J Specimen," 
by  L.  Nash  Gifford  (Sep  1978).  IC 
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EXPLANATION  OF  J (SINGULARITY)  RESULTS 


For  elastic  crack  problems,  the  asymptotic  expansion  of  the  elastic 
singular  solution  about  the  crack  tip  has  always  been  convergent  over  the 
entire  region  upon  which  it  has  been  imposed,  i.e.,  over  the  core  element 
or  over  the  enriched  element  The  higher  order  terms  neglected  in  the 
elastic  singular  solution  displacement  fields  become  increasingly  important 
with  distance  from  the  crack  tip;  they  can  apparently  be  accommodated  by 
the  regular  contributions  to  the  displacement  field  in  the  12-node  enriched 
element  case.  For  the  core  element,  on  the  other  hand,  their  importance 
has  been  minimized  by  keeping  the  element  radius  at  a very  small  value. 
Thus,  calculations  with  enriched  or  core  elements  in  the  elastic  case  have 
always  led  to  the  successful  prediction  (directly)  of  the  amplitude  of  the 
elastic  singular  solution  (K^  or  JSSy)- 

The  domain  of  convergence  of  the  plastic  singular  solution  in  elasto- 
plastic  crack  problems,  however,  is  limited  to  a region  well  within  the 
elastic-plastic  boundary,  i.e.,  where  the  plastic  strain  contributions 
dominate  the  corresponding  elastic  ones.  As  an  example  of  how  this  assump- 
tion is  easily  violated.  Figure  4 shows  the  near-tip  finite  element  mesh 
(the  region  within  1/2  in.  (12.7  mm)  from  the  crack  tip)  for  the  14-element 
enriched  plane  strain  test  case  of  Table  3a  with  n = 1 and  a = 10.  The 
enriched  elements  for  this  case  are  considerably  smaller  than  would  be  used 
in  the  elastic  case.  Superimposed  on  the  finite  element  grid  in  Figure  4 
are  the  elastic-plastic  boundaries  for  various  values  of  the  applied  remote 
stress.  It  can  be  seen  that  at  an  applied  remote  stress  of  20  ksi  (138 
Mpa),  the  plastic  zone  is  still  quite  small  and  occupies  an  area  less  than 
one-third  that  of  the  already  relatively  small  plastically  enriched  ele- 
ments surrounding  the  crack  tip.  At  higher  applied  stress,  the  situation 
improves  somewhat,  but  even  at  50  ksi  (345  Mpa),  one  of  tbe  enriched  ele- 
ments is  not  fully  within  the  plastic  zone,  and  neither  of  the  enriched 
elements  are  well  within  the  plastic  zone.  Recall  further  that  the  en- 
riched element  calculations  are  based  on  incremental  theory  plasticity. 

Thus  the  stress  state  at  the  higher  load  increments  is  dependent  on  the 
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Figure  4 - Elastic-Plastic  Interfaces  for  Fourteen-Element 
Enriched  Test  Case  of  Table  3 


22 
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stress  history,  i.e.,  inconsistencies  at  lower  amplitude  loading  may  prop- 
agate to  influence  the  50  ksi  (345  Mpa)  load  solution.  Were  the  ideali- 
zation one  of  plane  stress,  on  the  other  hand,  the  plastic  zones  would  be 
larger  at  a given  load,  but  still  not  sufficiently  large  to  enclose  the 
enriched  elements  well  within  the  plastic  zone.  The  natural  tendency  would 
be  to  overcome  this  difficulty  by  drastically  decreasing  the  size  of  the 
enriched  elements,  but  that  is  precisely  what  this  work  sought  to  avoid. 

g 

Moreover,  elastic  work  indicates  a loss  of  accuracy  as  enriched  element 
size  is  made  exceedingly  small. 

We  conclude  that  imposing  the  plastic  singular  solution.  Equations 
(2),  (9),  or  (10),  over  regions  which  violate  the  limitation  that  the 
plastic  singular  solution  be  contained  well  within  the  plastic  zone  leads 
to  incorrect  results.  Further  substantiation  of  this  conclusion  can  be 
found  by  considering  the  incremental  finite  element  calculations  involving 
enriched  elements  as  reported  in  Tables  2 and  3b.  At  low  values  of  ap- 
plied load,  the  plastic  zone  and,  therefore,  the  region  of  applicability 
of  the  plastic  singular  solution  is  small.  In  fact,  the  magnitude  of  the 
lowest  load  shown  in  Tables  2 and  3b  was  chosen  so  that  the  most  highly 
stressed  quadrature  point  in  the  enriched  elements  is  exactly  at  yield 
while  all  others  are  at  or  below  yield  stress.  This  condition  corresponds 
to  small  scale  yielding,  and  the  elastic  singular  solution  should  describe 
the  near  field  behavior  in  the  enricned  elements.  For  plane  stress  with 
V = 0.5,  the  asymptotic  elastoplastic  displacement  field.  Equation  (9),  is 
exactly  (1+a^)  times  the  elastic  solution.  (Similar,  though  more  com- 
plicated, relations  hold  for  V ^ 0.5  and  for  plane  strain.)  If  the 
elastic  solution  holds,  = Kj/d+a^)  and  J = (H-a^K^/E  = Jgs^/(l+aN). 

The  results  in  Tables  2 and  3b  at  the  lowest  load  (first  yield)  exhibit 
just  this  behavior  in  both  the  cases  of  plane  strain  and  plane  stress, 
i.e.,  Jg  -sr  / (1+a^) , and  thus  support  the  conclusion  that  the  enriched 
element  is  modeling  the  elastic  singularity  on  the  first  load  increment. 

As  the  load  increases  incrementally,  the  plastic  region  spreads  to  include 


the  entirety  of  the  enriched  elements  and  the  value  of  improves  relative 
to  Jp,  but  not  quickly  enough  to  give  accurate  values  within  the  practical 
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range  of  loading  considered  here.  Thus  the  inconsistency  between  assumed 
and  actual  asymptotic  displacement  fields  imposed  at  lower  load  levels  per- 
sists at  higher  load  levels  where  the  elastoplastic  displacement  field  is 
likely  to  be  more  correct. 

The  elastoplastic  calculations  involving  the  use  of  the  semicircular 

core  element  exhibit  less  discrepancy  between  J and  J values  than  the 

s p 

enriched  element  results.  This  observation  is  believed  to  be  a direct  con- 
sequence of  the  fact  that  the  core  elements  employed  were  significantly 

smaller  in  dimension  than  the  corresponding  enriched  elements.  With  this 

exception  noted,  the  poor  Jg  performance  of  the  core  element  is  explained 

in  like  manner. 


DISCUSSION  AND  CONCLUSIONS 

Calculations  to  obtain  values  of  the  J-integral  for  a sample  crack 
problem  have  been  carried  out  using  four  different  near  tip  models:  the 
core  element,  enriched  elements,  1/9-4/9  induced  singularity  elements,  and 
conventional  12-node  isoparametric  elements.  The  J values  were  determined 
in  all  cases  by  path  integration  (J^),  and  by  direct  calculation  of  the 
amplitude  of  the  crack  tip  singularity  (Jg)  for  core  and  enriched  elements. 
An  aim  of  this  effort  was  to  develop  a method  for  calculating  unique  J 
values  based  on  the  crack  tip  singularity  and  thus  avoid  the  necessity  of 
choosing  integration  paths  and  comparing  resulting  J estimates.  The  results 
presented  demonstrate  failure  in  reaching  that  goal;  rather,  the  most  re- 
liable J estimates  are  obtained  from  path  integration  based  on  finite 
element  solutions  employing  some  form  of  specialized  near-tip  elements. 

Further  verification  of  the  accuracy  of  path  values  of  J calculated 

using  enriched  crack  tip  elements  with  a multilinear  material  model  has 

been  obtained  by  comparison  with  J values  from  the  ASTM  round-robin 

16 


The  path  values  of  J^,  calculated  in  this  manner. 


plastic  crack  problem, 
were  in  the  central  region  of  the  band  of  results,  leading  to  increased 
confidence  in  the  use  of  the  technique.  As  expected,  the  corresponding 
singularity  values  of  J(Js),  as  calculated  with  enriched  12-node  elements, 
were  substantially  lower  than  the  path  values  and  were  outside  of  the 
range  of  other  reported  results. 
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In  retrospect,  the  fact  that  the  specialized  finite  element  calcula- 
tions give  accurate  path  values  for  J and  yet  are  unable  to  yield  accurate 
Jg  predictions  can  be  explained  as  follows:  The  path  independence  of  the 
J-integral  implies  a 1/r  singularity  of  the  strain  energy  density  function 
at  the  crack  tip.*  The  elastic-plastic  asymptotic  analysis  of  the  two- 
dimensional  crack  problem  for  either  a power  hardening  or  a multilinear 
material  model  discussed  and  referenced  earlier  yields  a plastic  singular 
solution  which  indeed  contains  a 1/r  singularity  in  strain  energy  density. 
The  region  over  which  this  singularity  dominates  is,  however,  dependent  on 
geometry,  material  behavior  (hardening),  and  applied  load  amplitude.  It 
is  this  dependence  (of  the  region  of  applicability  for  the  plastic  singu- 
lar solution)  that  has  led  to  the  difficulties,  discussed  earlier,  associ- 
ated with  imposing  the  plastic  singular  solution  on  a region  (element  or 
elements)  and  solving  for  its  amplitude  Jg.  On  the  other  hand,  the  J path 
integral,  which  is  a measure  of  the  amplitude  of  the  near-tip  fields,  is 
independent  of  these  size  restrictions  and  is  applicable  over  the  range 
from  small  scale  yielding  to  situations  involving  significant  plastic  de- 
formation. Therefore,  accurate  values  of  J CAN  be  obtained  by  using  the 
path  integral  approach  in  conjunction  with  specialized  crack  tip  elements 
which  impose  the  1/r  singularity  in  the  strain  energy  density  at  the  crack 
tip . 

Although  the  present  effort  failed  in  the  objective  of  calculating  the 
J-integral  directly  as  the  amplitude  of  the  plastic  singular  solution,  the 
use  of  plastic  singular  elements  at  the  crack  tip  was  found  to  be  bene- 
ficial for  the  evaluation  of  J about  paths  remote  from  the  crack  tip.  Thus 
the  only  real  limitations  of  the  present  effort  toward  the  prediction  of 
fracture  are  imposed  by  the  restrictions  on  the  use  of  the  J-integral  as 
a fracture  criterion.  Briefly,  the  J approach  is  limited  to  prediction 
of  fracture  Initiation  from  a preexisting  flaw.  The  J-integral  serves  to 
measure  the  amplitude  of  the  near-tip  field  and,  given  its  critical  value. 


*This  may  be  seen  by  considering  a limitingly  small  contour  or  path 
surrounding  the  crack  tip.  A lower  order  singularity  in  the  strain  energy 
density  would  lead  to  a limitingly  small  value  of  J on  this  path  and  a 
higher  order  singularity  would  result  in  increasingly  large  J path  values. 
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to  predict  crack  growth.  This  approach  clearly  presupposes  the  same  failure 
mechanism  in  all  situations  where  J is  applicable.  Thus,  the  scale  of  the 
plastic  deformation  prior  to  failure  must  be  limited  such  that  the  mecha- 
nism is  flat  fracture  and,  for  example,  excludes  ductile  tearing. 
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